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ABSTRACT 


The  marginal  success  of  the  several  high-resolution 
frequency-wavenumber  (f  -  k)  techniques  to  date  is  cited 
from  the  literature.  Their  ability  to  resolve  signals 
from  two  closely  spaced  sources  is  not  markedly  superior 
to  that  of  ordinary  beamforming.  Moreover,  such  non-linear 
techniques  yield  distorted  magnitudes  and  azimuths.  The 
ordinary  f  -  k  "spectrum"  is  shown  to  be  no  more  than  a 
1-signal  estimator,  and  the  existing  high  resolution  tech¬ 
niques  to  be  but  variations  of  that  1-signal  estimator. 

In  this  paper  the  notion  of  the  wavenumber  "spectrum"  is 
set  aside.  Instead,  by  analogy  to  the  1-signal  estimator 
(the  ordinary  f  -  k  "spectrum"),  a  linear  M-signal  estimator 
is  developed.  The  high  resolving  power  of  this  technique 
and  the  fidelity  of  its  estimates  is  demonstrated  theor¬ 
etically  and  by  computer  examples  both  real  and  synthetic. 


Introduction 


The  mathematical  development  of  linear  high-resolution 
frequency-wave  number  analysis  was  presented  in  the  Semi¬ 
annual  report  to  AFOSR  for  the  period  ending  1  March  1975. 
Since  that  time  the  software  implementation  of  the  theory 
has  been  accomplished  and  the  technique  has  been  applied 
to  explosion  and  earthquake  signals  recorded  at  the  Large 
Aperture  Seismic  Array  (LASA)  in  Montana.  The  application 
of  the  technique  to  real  data  is  the  basis  for  this  report, 
although  for  the  purposes  of  completeness  and  continuity 
the  previously-reported  section  on  mathematical  development 


is  also  included. 


The  frequency-wavenumber  spectrum,  which  is  a  multi¬ 
dimensional  equivalent  of  the  ordinary  frequency  spectrum, 
is  used  in  the  sciences  for  theoretical  and  experimental 
analysis  of  traveling  waves.  It  was  introduced  formally 
into  seismology  by  Burg  (5)  in  an  application  to  data 
analysis.  The  ordinary  unsmoothed  three  dimensional 
frequency-wavenumber  spectrum  of  time  series  data  sampled 
at  discrete  points  in  space  is  given  by 

M '  s  1  _  2- 

-jq-/C  {  AnCuJje^pEw^o,))]]  •  exp(ik-rj 

Cl) 

where 

A  is  the  index  of  the  spatial 

sample  points, 

AaM  explicxn(uj)] 

is  the  finite  Fourier  trans¬ 
form  of  the  nth  time  series. 

is  the  vector  wavenumber. 

is  the  vector  location  of  the  ntk 
sensor,  or  sample  point. 

Each  Fourier  transform  term  is  equivalent  to  a  sinusoid. 

For  example,  the  sinusoid  for  the  r\tk  transform  at  frequency 
Ol)  has  amplitude  Aa^^)  and  phase 
(at  the  center  of  the  time  window) . 

Now  —  k *  is  the  phase  delay,  between  the  origin 


k 

“n 


P(fe),k)  = 


1 


2 


and  r V\  ,  of  a  plane  wave  arriving  from  the  azimuth  of 

the  vector  k  and  traveling  at  the  phase  velocity 

V  - 

So,  multiplication  of  the  transform  by  the  kernel  eX.p  (l  k'l'^v) 
has  the  effect  of  advancing  the  sinusoid  by  just  the  amount 
the  wave  itself  had  delayed  it.  Thus  the  summation  in  (1), 
above,  is  a  beam  sum,  and  the  -P-k  spectrum  is  just  the 
frequency  domain  equivalent  of  ordinary  beam  steering. 

When  the  traveling-wave  delays  are  exactly  compensated 
for  by  the  beam  shifting,  i.e.,  when  the  true  k  of  the 
signal  is  selected,  the  sinusoids  add  up  constructively  with 
no  interference,  and  the  power,  P/  is  maximized.  Within 
certain  limits,  then,  maxima  or  peaks  in  P**k  space  are 
treated  as  indications  of  the  presence  of  traveling  plane 
waves,  and  the  location  and  size  of  the  maxima  are  taken 
as  estimates  of  the  speed,  bearing,  frequency  and  power 
of  those  signals.  If  more  than  one  signal  is  present  or 
if  there  is  noise  in  the  data,  though  exact  determinations 
are  no  longer  possible,  the  p-k  spectrum  is  still  useful 
for  detecting  and  estimating  signals,  again  within  limits. 

One  of  those  limitations  is  imposed  by  the  finite  width 
of  the  maxima  associated  with  signals.  (  9  )  The  case  is 
analogous  to  that  of  the  ordinary  frequency  spectrum  in 
which  components  are  represented  by  peaks  of  finite  width. 
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Plane  wave  signal  peaks  in  the  t-R,  spectrum  have  a  half¬ 
power  width  of  the  order 


_L 

AX 


where  AX  is  the  width,  or  aperture,  of  the  array  of 
spatial  sample  points.  (4)  If  two  signals  in  the  same 
time  window  and  frequency  band  are  also  close  enough  in 
phase  velocity  and  azimuth  so  their  wavenumbers,  say  k, 
and  kx  ,  are  such: 

|k,-kj  <  a  k  , 


then  their  maxima  in  the  f~K  spectrum  are  merged  and 
form  a  single  peak.  (23)  Thus,  because  the  sensor  arrays 
are  spatially  finite  their  resolving  power  is  finite.  Attempts 
to  increase  that  power  of  resolution  through  data  processing 
technique  have  required  mathematical  schemes  to  reduce  the 
width  of  the  lobe  of  the  signal  peak  (1-3,  6-15,  17,  19) . 

However,  the  straight-forward  geometric  appeal  of  this 
approach  has  proved  misleading  thus  far.  In  such  hybrid 
spectra  signal  lobe-widths  indeed  have  been  narrowed 
substantially.  Nevertheless,  when  signal  pairs  approach 
each  other  in  the  k-plane,  resolution  still  fails  as  the 
separation  nears  A  k,  to  wit,  the  lobe  half-width  for 
the  ordinary  -f-k  spectrum.  (2,  11,  13,  15,  20). 
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Observations  of  other  investigators  on  the  shortcomings  of 
various  high-resolution  frequency-wavenumber  techniques  are 
cited  below. 

Lintz  (15)  finds  that  the  high-resolution  f-k  spectral 
technique  of  Haney  (14)  does  not  significantly  improve  the 
capability  of  a  seismic  array  to  detect  multiple  time-overlap¬ 
ping  events  from  different  azimuths. 

Galat  and  Sax  (13)  experimentally  find  the  high-resol¬ 
ution  f-k  spectrum  of  Haney  (14),  and  that  of  Capon  (8),  (9), 

no  better  at  resolving  two  simultaneously  arriving  waves  than 
the  ordinary  f-k  spectrum.  McCowan  and  Lintz  (17)  call  atten¬ 
tion  to  an  unrecoverable  distortion  of  the  true  amplitude 
spectrum  in  Haney's  technique,  and  the  marked  disadvantage  of 
spurious  peaks  under  certain  conditions  which  they  regard  as 
the  inevitable  result  of  using  a  high-gain  procedure. 

Seligson  (20)  describes  conditions  under  which  Capon's 
high-resolution  technique  displays  less  "angular  resolution" 
than  ordinary  beamforming.  McDonough  (18)  concludes  that 
variations  in  amplitude  from  sensor  to  sensor  may  be  expected 
to  produce  anomolous  behavior  in  Capon's  processor.  Of  course, 
just  such  variation  in  amplitude  from  sensor  to  sensor  will 
result  precisely  because  of  the  presence  of  two  or  more  sig¬ 
nals.  McDonough  offers  arguments  to  show  that  ordinary  beam¬ 
forming  is  less  susceptible  to  instability  resulting  from 
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small  signal  modeling  errors  than  all  other  array  processors. 

Haney,  too,  notes  that  in  the  processor  he  describes 
(14)  variation  in  amplitude  from  sensor-to-sensor  could  dis¬ 
tort  the  spectrum  beyond  recognition.  He  remedies  this  diffi¬ 
culty  by  forcing  the  same  amplitude  upon  each  input  channel, 
thus  destroying  the  very  amplitude  information  that  would  be 
indicative  of  the  presence  of  two  or  more  signals. 

Woods  and  Lintz  (23)  conclude  that  given  favorable  con¬ 
ditions,  the  resolving  power  of  the  maximum-likelihood  f-k 
spectrum  can  be  effectively  infinite,  but,  disappointingly, 
offer  computer  examples  on  synthetic  data  in  which  the  input 
signal  pairs  are  well  spaced  to  begin  with  (they  are  separ¬ 
ated  by  a  distance  of  0.9  of  the  main-lobe  half-width).  Cox 
(11)  also  offers  theory  suggesting  that  given  arbitrarily 
high  signal-to-noise  ratios  arbitrarily  fine  resolution  should 
be  possible,  but  he  does  not  offer  a  method. 

It  may  be  argued  that  the  limited  resolving  power  of  the 
several  high-resolution  techniques  results  from  the  wavenumber 
spectrum  being  in  reality  a  1-signal  estimator.  Indeed  the 
ordinary  f-k  "spectrum"  is  a  least  squares  estimator  for  fit¬ 
ting  data  to  a  single  plane  wave,  as  shown  further  on.  In 
routine  automated  processing  of  the  LASA  LP  data  Mack  and 
Smart  (24)  found  the  ordinary  spectrum  useful  for  estimating 
only  one  signal  at  a  time.  Estimates  of  a  possible  second 
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signal  were  made  by  recomputing  the  wavenumber  spectrum 
after  the  first  (and  larger)  estimate  had  been  subtracted 
from  the  data.  They  call  this  process  stripping;  it  is  use¬ 
ful,  of  course,  only  for  estimating  signals  separated  by 
about  the  reciprocal  of  the  array  diameter  or  more.  At  that, 
such  estimates  of  a  pair  of  signals  are  not  optimum,  but  first 
order  approximations. 

Properly,  the  f-k  spectrum  is  defined  only  for  signals 
of  infinite  spatial  extent  traversing  infinitely  large  arrays. 
The  effect  of  a  signal  of  wavenumber  k  is  then  confined  to 
the  point  k  in  the  spectrum.  Approximations  to  this  defin¬ 
ition  are  useful  if  the  dimensions  of  signals  and  arrays  are 
sufficiently  large.  Failing  that,  the  "spectrum"  reduces  to 
a  1-signal  estimator  as  noted.  While  the  high-resolution 
techniques  do  attempt  to  extend  the  effective  array  diameter, 
they  all  test  the  wavenumber  space  with  a  1-signal  probe,  as 
in  the  ordinary  f-k  spectrum. 

It  is  proposed  here  to  set  aside  the  notion  of  a  spec¬ 
trum.  Rather  we  will  extend  the  1-signal  estimator  to  an 
M-signal  estimator  thus  to  permit  the  simultaneous  removal 
of  the  effects  of  one  signal  from  the  estimate  of  another 
and  so  achieve  true  high-resolution.  At  the  same  time,  use 
of  simple  beamforming  (in  the  k-plane)  to  estimate  each  of 
the  M  signals  will  preserve  the  stability  and  estimate  fidel- 


ity  of  the  ordinary  f-k  spectrum. 


In  the  following  discussion  a  1-signal  least  squares 
estimator  is  developed  and  is  identified  with  the  ordinary 
f-k  spectrum.  Analogy  to  the  1-signal  estimator  is  used  to 
develop  an  M-signal  estimator. 


Conventional  Frequency-Wavenumber  Analysis 


In  the  conventional  frequency-wavenumber  spectrum 
(ordinary  or  high-resolution)  a  single  plane  wave  is 
hypothesized  at  each  frequency.  That  model  is  then  tested 
over  the  wavenumber  space  of  interest.  One  attempts  to 
minimize  the  error 


<£~  =* 


IL  “A 


a. 


by  varying 


A  and  k 


where 


Un  are  the  complex  Fourier  series  terms 

(for  the  given  frequency) 

^  is  the  sensor,  or  channel,  index 

fsj  is  the  total  number  of  sensors 

are  the  location  vectors  of  the 
sensors 


A 

—  -W— 

k 


is  the  complex  Fourier  series  term 
for  the  hypothesized  plane  wave 
(at  the  given  frequency) 

is  the  wavenumber  of  the  hypothetical 
plane  wave  (at  that  given  frequency) 
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A/  is  the 


r\  =/,...  , 

model,  i.e.,  the  hypothesized  plane  w .we. 

Note  that  also  one  can  write  €  as 


a/ 

=  z: 

A*/ 


u 


-t  k  •  rk 


'  -A 


since 


=  / 


For  a  given  k,£  is  minimized  by  setting  A  to 

A  =  T^rZ? 


a**/ 


which  is  shown  by  the  following: 

c9„  +  i  cu  s  Ua  e"'^^ 
<3  -tic.  S  A 


Let 

and 


Then 


A/ 


A  =  / 

aJ 


(aA-e) 


=  z  (^  -  ccA-cr 


A®  t 


Take  partial  derivatives: 

A/ 
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<3 


■  Sr#*  =  ZA/ 


Setting 


=  <P_£ 
«?C 

A/ 


_  /  y' 

A/  4^ 


C  * 


A  =  / 


/  44 

ttz:c- 

«=/ 


and 
<£>  -/ 


;c  =  4  =  =  77 4  tJ„e' 

'  ^  r\a/  A* / 


So,  minimized  with  respect  to  A, 

*  =ZT  |  -  ij  Z  U7- 

7  /»/  J 


A-/ 
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This  expression  can  be  separated  into  2  parts,  thus: 


^  ~  (c5^  —  <3  j  4  ("C  ~ 


As/ 


a7 


A/ 


As  / 


Thus, 


A/ 

z; 

A  «  / 


a/ 


A  =■  / 


c3  </■  c3  * 

+ 

-  £c2  -c5  A/  c9 

. 

z 

t  CA 

z 

1 

uA 

~  M 

A/ 
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A=  / 
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The  second  term  is 
spectrum 


the  ordinary 


P(f,k)  = 


/ 

A/ 


frequency -wavenumber 

A/  7*— 

1.  Ik 

=  /  J 


Since ^  is  a  squared  modulus 

>  o 


and 


since  it  is  a  sum  of  squared  modulii. 


S imilarly 
Since 


Pc  h  ^  ° 


A-/ 


A  =  / 


Z3  UJ  -P<b>oy 


a/ 


ziuT» 


So  to  minimize  ^  one  must  maximize  Pth 


a/ 


r-0  -  £  U*  -  A 

/W 
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becomes  exactly  zero  when 


U,  =  •  77  U,e'^5  =  Ae'^'5 


t\  =  3  y  ...  y  aJ 


that  is,  when  the  data  describe  a  single  plane  wave  exactly. 


The  smaller 


^  is,  in  a  given  situation,  the 


more  likely  is  the  hypothetical  plane  wave 

A  e-‘  k'r 

because  the  smaller  <£rC  k)  is,  the  larger  the  F-statistic 
is  for  the  hypothesis.  The  F-statistic  is  given  by 


F  m  ( M -i) •  PcT^ct) 


This  single  plane  wave  model  is  often  applied  in 
attempts  to  analyze  a  2-signal  case  (or  a  possible  2-signal 
case) .  In  such  an  analysis  each  signal  is  treated  as  if  it 
existed  by  itself,  the  presence  of  the  other  being  ignored 
with  consequent  distortion  of  estimates  by  mutual  interfer¬ 
ence.  This  interference  can  be  serious,  and  if  the  two 


signals  are  not  separated  in 


-space  by  at  least  the  half¬ 


width  of  the  main  lobe  of  the  array  response,  they  are  likely 
to  appear  as  but  one  signal,  their  main  lobes  having  coalesced. 
Attempts  to  improve  the  performance  of  the  single  wave  hypo¬ 
thesis  (in  application  to  the  two  signal  case)  have  been 
made  in  which  the  main  lobe  of  the  array  response  has  been 


slenderized  mathematically  by  alternative  methods  of  esti- 


mation  of  the  wavenumber  spectrum.  The  object  has  been  to 
reduce  the  main-lobe  half-width  and  so  resolve  signal  pairs 
which  otherwise  have  coalesced  main-lobes  indistinguishable 
from  a  signal  case.  These  results  have  been  marginal.  In 
the  various  high-resolution  techniques  the  influence  of  the 
one  signal  on  the  analysis  of  the  other  has  been  ignored. 

Analysis  of  possible  2-signal  cases  calls  for  a 
2-signal  model,  in  particular  when  the  2-signals  are  known 
(or  suspected)  to  be  so  close  together  as  to  have  their  main 
lobes  merged. 

As  the  1-signal  model  serves  for  both  the  0-  and  the 
1-signal  case,  so  one  might  expect  a  2-signal  model  to  be 
effective  in  all  three  cases:  0,  1,  or  2-signals. 


Multiple  Signal  Frequency-Wavenumber  Analysis 


By  analogy  to  the  1-signal  model,  one  would 
expect  to  solve  a  2-signal  model  by  minimizing  the  error 


A/ 


*  -  z: 


n  =  / 


UA-  Ae,^''s 


varying  A,  k,  B,  and  K<s  ,  where 


B 


k 


is  the  complex  Fourier  series 
term  for  the  second  hypothesized 
plane  wave  (at  the  same  given 
frequency) 

is  the  wavenumber  of  the  hypothe¬ 
tical  plane  wave  (at  that  same 
given  frequency) 


There  are  now  two  signals  to  solve  for: 


A  e 


i  k  •  r„ 


r\  d  E> 


Let 


"T;  =  LL  -  A  e 


1  k  *  Ky 


-  E> 


i  k 


then 


Again,  let 


=  z 

* -f 


x 


/V 


=  Z  T/T, 


n  =  / 


A  *  *  /  c  ,  A 


c a  -  /  c 


Taking  first  partial  derivatives  while  noting  that 

£  A _ 


a  A  a  A*  _  i 


and 
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ac 


a  c 


/  > 


<3 

<P  <3 

A/ 

-  £  T*(- 

r\  */ 

and 

A/ 

<?<£ 

PC 

-r  t*(- 

n  =  1 

Setting 

36  _ 

~  P  C 
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t  Ic-’r' 


)  +-r;(-e-iUv~) 

ik'^)+  T(erlUt) 


as  in  the  1-signal  case, 


^  ^  ^  i  <5 


=  O 


A/ 


■  ■«  «  Tv 

Therefore, 


^  = 


>n=r/ 


A  =  ^  (llA  -  E>  ^'r* 

*\ = /  / 


Analogously 


*\=! 
A/ 


&  -  A/  ~~  A  ^  ^  ’  er~ 1  ^  ■  r- 


In  this  form  A  and  B  are  optimized,  that  is,  they  produce 
the  minimum  value  of  ^  for  any  arbitrary  pair  of  and 

/i*.  Adopting  the  notation: 

p-ttZ:  Uv^t 

M 


Q*  n  22  ^ 


Ae  I 

N 


E -TirC  *' 

l\c| 


ick-IcK 


one  may  write  simply: 


16 


A-P-B-E  **=1  B=  Q-A-E* 

Rearranging  to  solve  A  and  B  simultaneously 

P  =  A  +  BE 
Q=A-E,‘-t  E> 


p 

E 

1 

P 

a 

1 

_  ,  R  -  . 

E* 

a 

i 

E 

1  D 

1 

E 

E* 

1 

E* 

1 

A  =  (P-Q  E)/(  t  -  E*E) 

B  =  (Q-PE*)/(|-E*E) 


Written  out  at  length, 


and  e>  is  similar  in  form. 

Introducing  a  factor  of  y.  into 

Nl 


nr  t 


A*  I 


N 


“  kT  c  (  u  *  -  A*e-'  '‘••"A  -  B*V'"k  •>''.) 

'  .iXe 


x(ll„-A 


-B 


m 


-(  B*tt  t  +  B  i  £  u;  -1 

'  ^  =  1  r\tl 

+  (  A* A  •+  B*  B) 

+(a^Bttz:  e,(1-X)F"  +ab^tVz:  •-*  cl-lR) 

*=  l  l\S| 

*r  =  -itZ:  U*U*  -  (  A*P  +  AP*)  -  (B*C)l  +  BQ*) 

^  *N=i 

•\ 

+  (A*A4B*B)  +(a‘BE  +  AB^E") 

Rearranging  the  terms  in  ^ , 

=  -t  Z  u*u„  -  (  a*  p  +  ap*)  -  (  B*Q  +  B  a*1) 

™  A=l 

+  A*(A4BE)  +b*(ae*  +  b) 

and  recalling  that 

p  =  A  +  BE  and  Q -  ae*  +  b 

«  -  Triu^u,  -(AP*  +  BQ*) 
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Further,  substituting 

A  -  P-  BE  and  B~Q.-AE* 

-  -jj-il U-J* -(pp+q* o  -  PQ*E*-P*QE) 

*'»•  .  l—  H-  r- 


or,  written  out, 


The  identity  of  these  last  2  equations  may  be  demonstrated 
by  noting  that  the  numerator  (above)  equals 


.\=i 


=  P*P+Q*Q  -PQ*E*-P’,QE 

Since  ^  is  a  sum  of  squares,  by  definition  it  must 
be  non-negative  everywhere.  Therefore  the  second  of  the 
2  terms  in  ^  ,  above,  must  always  be 
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^  n  2juj 

.N  -I  1 

Thus  to  minimize  ,  one  must  maximize 


»\  =  i 


This  is  the  2-signal  test,  analogous  to  the  ordinary  frequency- 
wavenumber  spectrum,  which  is  the  1-signal  test.  However,  it 
is  more  convenient  to  retain  the  form 

AP*  +  BCT 

This  2-signal  -f"“  k  "spectrum"  then  is  computed  from  3 
beams  (as  the  ordinary  spectrum  is  computed  from  1  beam) . 

The  beams  are 

the  mean  of  the  data  transforms 
that  have  been  beamed  to  k 
(one  of  the  two  wavenumber 
variables ) , 

the  mean  of  the  data  transforms 
after  beaming  to  (the 

other  wavenumber  variable) , 


which  is  the  (complex)  array 
response 

This  2-signal  test  is  solved  as  is  the  ordinary  -f-k 


spectrum,  numerically,  by  searching  the  wavenumber  space  of 


•do 

interest.  Now,  however,  there  are  4  dimensions  to  search, 
over  which  to  test  the  error  criterion. 

It  is  instructive  to  submit  a  known  pair  of  pure, 
noiseless  signals  to  the  2-signal  test  to  illustrate  the 
function  of  the  elements  of  the  expression: 


U  -  F 


inTiH  SSiS^,?’  k  and  *  ,since  these  are  kno“" 


P  =  TxT <G  (  F e’  k'r'v  +  G  <sl 

-\=l  1 

~  nT  C  P  +  <G 


.-i  k: 


INI 

F  +  Q~bZ* 


.Ok-k) 


=  F  +  GE 


ct7'  .  r;  iFik\  ,,-i'kp 


Q  =  wZlCFeP-^  +  q 

r\=  |  1 

-  FE*  +  q 


)•  e= 


Then 


(P  +  qe)/o-e^) 

(f+qe  -cfe*+q)e)/(i-e*e) 

(  F  +  QE  -  FE^e-Q£)/(i-e^) 

F(|  ~  E*e)/(i~e:*e)  =  F 


B~(h  L4  +  q  -(F  +  GE)E*)/(|  -E*E) 
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B  =  GO  -E^E^I-E^e)  =  G 


Thus 


ki 


=  ttz: 


►\=i 


U»  -  A. 


i  Lr« 


"  -ge'U 


nJ 


t  k- 


4Qe' 


iH: 


F< 


i  krk 


&  S  O 


This  little  exercise  clarifies  a  bit  the  function  of  the  array 
response,  E,  in  the  signal  models  A  and  B 

The  development  of  the  2-signal  test,  of  course, 
suggests  the  derivation  of  a  3-signal  test,  by  analogy: 

First,  the  form  of  the  test  would  be,  analogously, 


Ni 

I 

N  -O 

r\- 1 


U«  -  A 


k- 


Wv 


-B 


i  k: 


c 


it- 


Introducing  the  notation 

R  ^  -hZ  U* 

A=l 

r- _ L  i  (‘fe  -  k)’^ 

As  I 

and  expanding  <£T  : 


e.^z: 


N  itk-k)-^ 


A=1 

NJ 


Ej  =  -ntZZ  =' 


iC^-kVr* 
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-( u*e-iU  +  a  -kz  u;^ k~-) 


■(B*^u^  +  BAZu:^*) 

^ "  I  r\c|  ^ 

+  c^2u^-7') 

+  ( A*A  +  B*  B  +  c*c) 


r\=t 


4  (  A'Bxri  e,(4-kK  +  A  BMri  ei(t'  '‘H 

**•  A*.  ' 

+  ( B’CA-Z  B  C^Z>i(W*) 

^  *\-(  7 


•V 

M 


4  (  A*C-srZ  el(t'L>^  +  A  C%t 

A=l  «\-l  1 


M 


=  n2  u,*u 


rv=  i 


A  W/y 


-  (  A*  P  4  RP‘)  -  (B*Q  4  BQ*)-(c*R4  CR‘) 
4  (  A*  A  4  B*B  4  C*c) 


4(A*BE,+AB*E,*)  +  (B*CE3+  BC*Ef)  *  (A*CE,+AC*E0 


Now  noting  that  in  the  2-signal  test 


P  =  A  ■+■  B  E  and 

Q  =  AE*+  e> 


so  that 


A 


P  E 
Q  1 
I  E 
E*  I 


I  P 
E*  Q 

I  E 
E*  I 


one  recognizes  that,  in  the  3-signal  test. 


P  =  A  4  B  E,  4  CE» 
Q  =  AE,*+  B  +  CE3 
R  =  A  E*K  +  B  Ej  +  C 


and,  defining 


cJe  ^ 


I  E,  E* 
E*  I  E . 
Ez*  E*  I 


A  -  [  P(l  -E‘E3)  4  Q(e*E,-E,)  4  R(E,Es-E^]/clen 
B  -  [  P(E*E4-ED  +  Q(|-Ex*Ej  4  R(E^ErE^cJer, 
c  -  [  p  ceTe^-et)  +  Q(e;e1*-e;)4  r(i-e,*e)]A 

d«s.fs=  1-EfE,  -  E*E»  -  EJEj  4E,E*ES  +  E*  E,  E* 
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Now  rearranging 


I  N 

nZ  uX-  (a*P+  b*q+c«r)  -  (ap*  +  Ba^  +  CR*) 

•\«<  ' 

4-  A*(A+  BE,  4  CE,)  +  ^(AEf+  B  +  C  E3)  •+C<*(AE^: 4  BE.*  +  C) 

and  substituting  P  ,Q  ,  and  R 

I  , 

<f=  nZ  UJlh  -(AP*4  BQ*  4  CR*) 


To  minimize 


,  then,  one  will  maximize 


AP“  +60*+  CR* 


the  3-signal  test,  or  3-signal  analog  to  the  conventional, 
1-signal  frequency-wavenumber  spectrum.  The  function  is 
composed  of  6  beams:  P  .Q  ,  and  R  ,  the  3  beams  of  the 
data,  U,  ,  and  EE,  ,  EE* ,  and  EE3 ,  the  3  beams  of  the  array 


response . 


Remembering  the  1-signal  test  (conventiona 


i  f-k 


spectrum) , 


-  tZ  U»-  A 


njZ)  (u  *-  A*  e'*  -Ae‘  ^ 


we  may  rewrite  it  as 

i  i  4^ . ,  « 


=  nZ-  u*  u«  -  a*p  -  a  p*  -  a* a 
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(since  Ps  N  ZI  ^  **  )  . 


or 


w 


«  =  7J 72  UfU.  -  A*A  -  A  P*  +  A*A 


A*l 

M 


=  wz:  uru,  -ap 


As) 


Thus 


AP*  ■  i(i|ZU,.:t  T) 


is  the  expression  one  must  maximize  in  order  to  minimize 
the  error.  So  the  "f-k  spectrum  (for  the  1-signal,  conven¬ 
tional,  case)  is 

AP*  , 

AP*  +  BCl*  is  the  2-signal  test, 


and 


and 


AP*+  BQ*  +  C  R*  is  the  3-signal  test 


In  the  1-signal  test 


a=_p 

M  I 


For  the  2-signal  test 

IP  E,| 
a  i 


i  e, 

Ef  I 


For  the  3-signal  test 


A 


B  = 


P  E,  E, 
Q  I  E, 
R  E?  I 


I  P 

E*  Q 


IE, 

Ef  I 


I  EE 
ETI  E 
EtEr  i 


1 


e.lc 


I 
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This  formalism  makes  evident  the  relationship  between  the 
successive  tests.  Thus  one  may  extropolate  and  directly 
write  the  expression  for  the  M-signal  test  in  simple,  terse 
form.  For  example,  the  4-signal  test  is 


Ap%  BQ*4  CR*+DS* 


in  which 


S, 


the  sum  of  the  data  beamed  to  yet  a  4th  point 
,  is  introduced  into  the  sequence  P.Q  ,  and  R  ; 


and  in  which 


P  E, E»E„ 
Q  I  E5E5 

r  e;i  e, 
s  e;e:  i 


I  e.e.e, 

E,  I  E5E5 

e;e:  i  e 
e;e;e:i  , 


and 


E<  is  the  array  response  at  (*-to .  E  g  ,  that  at 


(  —  k)  ,  etc. 

Note  that  the  four-signal  test  is  computed  from  10 
beams;  4  beams  of  the  input,  U*\  ,  and  6  of  the  array  response. 
In  general,  the  M-s  ignal  test  requires  M  beams  of  input 
data  (U)  ,  and  m(m-i)/2  beams  on  the  array  response,  for  a 
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total  of  m(m + \)/£  beams  to  compute  the  least-squares  error 
at  any  point  in  the  2M~dimensional  space.  But  the  beams  on 
the  array  response  are  computed  from  the  same  complex  trig¬ 
onometric  terms  that  are  required  for  the  M  beams  of  the  input 
data.  So  the  M-signal  test  requires  evaluation  of  z  MN  sine 
and  cosine  terms  to  compute  the  error  at  any  point  (  NJ  is  the 
number  of  sensors  in  the  array) .  Thus  the  number  of  trigo¬ 
nometric  terms  requiring  computation  increases  linearly  with 

M. 

It  must  be  noted  that  a  multiple  signal  test  is  not 
everywhere  well-behaved,  but  has  a  singularity.  For  example, 
in  the  case  of  the  2-signal  test,  if 


so  that 


and 


k 

Q— P 

E— 1  _ 

dzz-  is  undefined.  The  value  it  will  take  on  at  —  1^ 

depends  on  the  direction  from  which  fe—U.  Though  this 
can,  of  course,  be  shown  analytically,  it  is  a  bit  tedious  for 
repetition  here.  The  contoured  map  of  an  example  (figure  1  ) 

displays  this  characteristic  graphically.  The  contoured 
function  is  the  2-signal  test 


A  P  +  BQ' 


-ffcr* 


with 


held  fixed  as 
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1^  held  fixed  as  *i<  varies  over  the  plane.  Note 
that  the  contour  lines  all  run  together  at  as  |<^ 

may  range  arbitrarily  close  to-  l<.  but  must  not 
take  on  that  value  exactly.  The  data  in  this  figure  consist 
of  2  closely  spaced  signals.  The  fixed  vector,  1*.  ,  was 

set  at  the  peak  of  their  merged  main  lobes. 

One  might  dismiss  this  singularity  from  practical  con¬ 
sideration  since  signals  of  identical  speed  and  bearing  are 
indistinguishable  by  array  methods.  The  test  for  2  signals  at 
the  same  wavenumber  location  is  thus  unnecessary  anyway.  But 
if  the  2-signal  test,  say,  is  applied  to  data  composed  of  only 
1  signal,  must  not  both  the  probe  vectors  approach  the  same 
point,  i.e.,  the  wavenumber  location  of  the  input  signal,  in 
order  to  merge  and  reduce  the  function  to  the  1-signal  test? 

We  have  seen  that  when  the  data,  ,  consist  of  the  same 

number  of  signals  as  that  for  which  one  is  testing,  the 
test  performs  as  expected:  the  error  _is  minimized  at  the 
wavenumber  location  of  those  input  signals,  and  the  signals 
are  recovered  undistorted.  Suppose,  though,  that  the  2-signal 
test,  say,  is  applied  to  data  consisting  of  just  plane  wave. 
Let 


U.  =  F 


i  It  •  rk 


TA  -  U,-Aeikrr'  -  be' 


ifc- 


in  the  error  expression 


-  Trim* 

*.*-1 


We  have  to  maximize 


AP%  BQ* 
k)*T^ 


P  =  F  TT  z:  = V ^  «J  Q  -  F  -kZ 

■  _  l\st  _ k  A=l 

If  k.  goes  to  ^  ,  then 

P  -  F,  Q  -  FE* 


and 


A  =  (P-QE)/(l-E*E) 

becomes 

A -(F-FE"E)/(l-E*E)  -  F 

and 

B  =  (Q-PE*>/(I-E*E)  =  ( F  E*-  FE*)/(l-E*E-) 

-  O 

and 

*~-bi  uru.  -  (AP*+  BQ*) 

-  (F'e^HFe^-*)  -  (F“F  +  O) 

=  iiz  F*F  -  F*F  =  o 

When  U.  goes  to  the  error  is  minimized,  the 

signal,  F"  ,  is  recovered  undistorted,  and  the  hypothesized 


second  signal  vanishes.  This  solution  is  invariant  though 


k. -plane. 


/l<  be  permitted  to  range  over  the  entire 
excepting  the  point  U  .  Thus  the  2-signal  test  does 

r  i 

spectrum  in  the  presence 


not  reduce  to  the  ordinary  "f  ”  k, 
of  a  single  plane  wave,  and  'fe.  is  not  required  to  go  to 
l<  nor  would  the  gradient  of  ^  with  respect  to  /tk. 
lead  to  l<  (if  one  were  using  a  steepest  descent  technique 


to  minimize  £  )  . 


Numerical  Solution  of  the  Multiple  Signal  Test 

One  might  propose  to  carry  out  the  numerical  solution 
of  a  multiple  signal  test  by  a  straightforward  search  of 
the  entire  wavenumber  space  of  interest,  as  is  done  in  the 
computation  of  the  conventional  f-k  spectrum.  But  the 
multiple  signal  test  may  be  used  in  more  practical  fashion, 
with  greater  efficiency,  as  a  follow-up  to  the  ordinary 
f-k  spectrum.  Since  a  high-resolution  array  process 
by  design  is  intended  to  separate  signals  otherwise  unre- 
solvable,  there  is  sound  justification  to  limit  its  use  to 
the  vicinity  of  signals  tentatively  identified  beforehand 
by  less  powerful  but  faster  techniques.  This  is  an  advan¬ 
tageous  circumstance,  since  an  M-signal  test  is  a  function 
of  2M  dimensions  of  wavenumber  and  would  otherwise  prove 
computationally  less  efficient.  Applying  the  2  -  signal 
test  to  the  highest  peak  of  an  ordinary  f-k  spectrum, 
then,  one  hypothesizes  the  presence  of  2  plane  waves  which 
appear  as  only  1  because  of  their  proximity.  By  the  hypothe¬ 
sis  the  spectral  peak  lies  within  the  area  of  the  main  lobe 
of  either  signal  and  thus  £  may  be  minimized  directly  by 
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3  2 


the  method  of  steepest  descent.  This  is  the  procedure  used 


here . 


Since,  as  has  been  shown  earlier, 


/k,  =  k 


is  prohibited,  the  descent  cannot  begin  from  any  one  single 
point  in  the  k  -plane,  as,  for  example,  the  peak  under  con¬ 
sideration.  But  any  pair  of  points  in  that  vicinity  is 
suitable;  all  lead  to  the  same  solution.  A  convenient  pair 
are  (1)  the  peak,  and  (2)  the  adjacent  minimum  of  with 


respect  to,  say. 


when  K  is  fixed  at  the  peak  as  in  the 


previously  discussed  figure  1  .  The  gradient  of  £  is  com¬ 
puted  at  this  pair  and  ^  itself  then  recomputed  at  a  new 
location  down  the  gradient.  The  length  of  this  first  step 
in  the  descent  is  some  fraction  of  the  width  of  the  array- 
response  main-lobe,  thus  chosen  to  ensure  that  the  process 
does  not  jump  from  the  vicinity  of  the  solution  into  the 
range  of  an  adjacent  relative  minimum.  The  gradient  is  new¬ 
ly  computed  at  this  second  location;  another  somewhat  smaller 
step  is  taken  down  the  gradient;  the  gradient  is  once  more 
computed,  now  at  this  third  location,  and  so  forth  in  suc¬ 
cessively  smaller  steps  until  the  point  is  reached  in  that 


4-dimensional  space  at  which  the  gradient  goes  to  zero. 
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Some  examples  with  synthetic  data  for  the  LASA  LP  array 
follow. 

Figure  2  represents  the  array  response  of  LASA.  The 
function  is  mapped  in  contours  of  3  db  intervals  down  from 
the  peak  at  the  center.  At  present  7  long-period  vertical 
seismometers  remain  at  LASA:  the  D-ring,  2  elements  of  the 
C-ring,  and  the  center  sensor  at  AO.  Thus,  the  half-width 
of  the  main  lobe  is  about  0.016  cycle/km.  At  20  seconds 
period  and  3.5  km/sec  that  half-width  intercepts  nearly  70 
degrees  of  azimuth. 

Since  the  error  expression  for  the  2-signal  test  is  a 
function  of  4  dimensions  of  wavenumber  and  cannot  be  presented 
in  map  form  as  are  ordinary  f-k  spectral  sections,  numerical 
presentations  must  be  resorted  to.  Figure  3  presents  the 
first  and  last  page  of  computer  print-out  of  the  successive 
steps  in  the  solution  of  a  single-frequency  test  case  cf  2 
noiseless  signals  separated  by  1/8  of  the  array  main-lobe 
half-width.  One  signal  is  2  magnitudes  greater  than  the 
other  (ten  times  the  amplitude)  and  is  180°  out  of  phase 
with  it  at  the  array  center  such  that  they  destructively  in¬ 
terfere  with  each  other.  At  the  upper  right  of  the  first  page 
are  the  signal  descriptions;  beneath  that  are  the  array  coor¬ 
dinates  (in  km.)  and  the  Fourier  transforms  of  the  signals. 
Thereafter  each  successive  line  represents  a  step  in  the 


m 


tit 
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optimization.  The  E  format  number  at  the  left  in  each  line 
is  the  error  term.  The  amplitude,  phase  and  coordinates 
(kx,  ky)  at  each  pair  of  points  are  given  in  the  two  columns 
enclosed  by  vertical  lines.  The  unit  vector  of  the  gradient 
is  given  by  the  4  numbers  at  the  right  of  the  page.  The  size 
of  the  step  from  the  previous  point  just  precedes  the  unit 
vector . 

The  solution  is  given  at  the  bottom  of  the  second  page 
by  the  complex  transforms,  wavenumber  locations  and  the  final 
error  value.  The  solution  is  both  accurate  and  precise;  the 
high  resolution  has  introduced  no  distortion  such  as  charac¬ 
terizes  the  non-linear  techniques. 

The  computer  routines  of  figure  3  that  apply  the  2- 
signal  test  were  introduced  into  a  general  f-k  analysis  pro¬ 
gram  called  FKSCAN  which  was  styled  after  FKCOMB  [jVlack  and 
Smart  ( 24 )j  .  To  this  automated  high-resolution  processor 
synthetic  time  series  were  submitted  for  trial  analysis.  One 
test  consisted  of  a  unit  plane-wave  from  356°  at  3.5  km/s  to 
which  synthetic  random  noise  was  added  to  make  the  signal-to- 
noise  power  ratio  equal  to  4 .  To  this  combination  was  added 
a  second  plane-wave,  2  magnitudes  larger  than  the  first,  from 
302°,  also  at  3.5  km/s.  In  the  band  of  interest  '  5-23  seconds 
period)  the  2  signals  overlap  in  wavenumber  space.  At  23.3 
seconds  period  they  are  separated  by  0.7  the  main-lobe  half- 
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width;  at  16.0  seconds,  by  1.0.  Each  signal  was  of  20  seconds 
period  enclosed  in  192  second  cosine  envelope.  They  arrive 
at  the  array  center  at  the  same  instant.  A  256-second  time 
window  was  applied  for  the  analysis.  The  resulting  bulletin 
from  FKSCAN  is  given  in  figure  4,  and  is  self-explanatory. 

The  larger  signal  is  shown  arriving  from  301°  at  3.518  km/s ; 
the  smaller  from  354°  at  3.059  km/s.  They  differ  in  apparent 
magnitude  by  2.05  (from  the  ratio  of  the  power  summed  over 
the  band) . 

The  last  item  on  the  second  page  gives  a  measure  of  the 
assurance  one  would  have  had  of  the  validity  of  such  a  de¬ 
tection  had  it  appeared  in  processing  of  real  data.  In  routine 
processing  of  such  4  minute,  16  second  windows  each  interval 
yields  2  suits  of  vectors,  or  detections.  At  the  rate  of  2 
suits  of  random  vectors  per  time  interval,  so  anamolous  an 
angular  concentration  of  vectors  would  appear  only  once  every 
11  days,  on  the  average.  |^This  algorithm,  installed  in  FKSCAN 
to  provide  an  independent  detection  statistic  separate  from 
the  F-statistic,  is  based  on  a  probability  expression  developed 
by  the  author  which  he  intends  subsequently  to  submit  for 
publicationTj 

Thus,  in  this  modest  test  of  the  2-signal  detector 
functioning  in  the  presence  of  random  noise,  the  small,  "hidden" 
signal  is  recovered  as  a  strong  detection. 


The  array,  the  relative  magnitudes  and  azimuthal  spac¬ 
ing  of  these  2  test  signals,  and  the  frequency  band  in  which 
the  search  was  made  anticipate  a  test  on  real  LASA  data  in 
which  similar  conditions  were  expected.  The  2  seismic  events 
sought  in  the  real  data  are  recorded  in  the  U.S.  Department 
of  the  Interior's  Earthquake  Data  Reports  36-74  and  43-74. 

On  31  May  74  at  0313:11  an  earthquake  occurred  in  the  vicinity 
of  Unimak  Island  in  the  Aleutians.  At  0326:57  a  second  event 
occurred  in  eastern  Kazakh  SSR.  The  first  quake  had  a  body 
wave  magnitude  of  4.8  and  a  surface  wave  magnitude  of  4.6. 

The  second  event  was  recorded  as  Mfi  5.9,  with  measurements  un¬ 
available.  The  Unimak  signal  was  expected  at  LASA  about  0330 
from  302°,  with  the  Kazakh  signal  expected  about  0406,  in  the 
ongoing  coda,  from  356°.  Figure  5  displays  the  seismograms 
of  this  interval  for  all  seven  LASA  stations.  The  figure  begins 
at  0329:08  and  continues  past  0415.  The  anticipated  onset 
of  the  Kazakh  surface  wave  is  marked  by  the  arrow  at  0406:09. 
The  circled  numerals  at  the  bottom  of  the  figure  number  the 
successive  time  windows,  indicated  by  arrows,  that  were  sub¬ 
mitted  to  the  high-resolution  array  processor. 

Figures  6  through  11  are  the  resulting  bulletins  for 
the  6  time-windows  marked  on  the  seismograms.  The  relatively 
narrow  band  from  16-23  seconds  period  was  chosen  for  this 
analysis  because  it  was  anticipated  that  the  faint  signal  from 


from  Kazakh  was  most  likely  to  appear  in  these  frequencies 
if  at  all.  As  before  there  are  two  sets  of  tentative  detections 
in  each  time-window,  that  is,  2  detections  per  frequency.  The 
suit  at  the  upper  left  contains,  at  each  frequency,  the  sig¬ 
nal  pick  of  greatest  power.  One  might  call  these  the  primary 
detections.  The  suit  at  the  upper  right  contains  the  smaller 
signal  picks.  Before  each  suit  is  submitted  to  azimuthal  dis¬ 
tribution  analysis,  the  program  computes  the  straight  line 
through  the  frequency-wavenumber  origin  (0,0,0)  which,  in  the 
least  squares  sense,  best  fits  all  the  vectors  in  the  suit. 

The  back  azimuth  and  phase  velocity  of  that  mean  are  then 
printed  in  the  bulletin. 

The  back-azimuth  of  the  mean  of  the  primary  signal  picks 
in  succeeding  time-windows,  then,  for  these  6  intervals  reads, 
in  sequence:  328°,  319°,  323°,  353°,  356°,  and  318°.  In  win¬ 
dows  1,  2,  and  3  the  detector  is  "triggering"  on  the  ongoing 
coda  from  the  Unimak  earthquake.  But  in  windows  4  and  5  it 
turns  and  indicates  the  back  azimuth  of  the  Kazakh  site.  Then 
in  window  6  it  returns  to  the  Unimak  coda. 

Figures  12  and  17  (one  for  each  of  these  6  time-windows) 
are  contoured  plots  in  3  db  intervals,  of  the  conventional 
wavenumber  spectra,  integrated  over  the  frequency  band  0.043- 
0.063  Hz  (16-23  seconds  period)  after  the  secondary  detections 
derived  from  the  high-resolution  processor  have  been  filtered 
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out.  These  spectra  make  visible  the  observations  of  the 
previous  paragraph. 


Conclusions 


The  high  resolving  power  of  the  linear  multiple-signal 
analysis  and  the  fidelity  of  its  estimates  have  been  de¬ 
monstrated  by  computer  examples  and  by  application  to  real 
signals . 

Computer  examples  indicate  that  this  technique  is 
capable,  in  the  absence  of  noise,  of  exactly  recovering 
the  amplitude,  phase,  and  velocity  of  two  simultaneously 
arriving  Rayleigh  waves  at,  for  example,  LASA,  which  differ 
in  azimuth  by  as  little  as  8°,  even  if  one  signal  is  10 
times  larger  than  the  other.  In  the  case  of  the  simultan¬ 
eous  arrival  of  a  small  signal  with  S/N  of  2  and  a  signal 
100  times  larger,  with  a  difference  in  azimuth  between  the 
two  of  54°,  the  magnitude  (Ms)  of  the  small  signal  can  be 
recovered  with  less  than  3%  distortion. 

The  extraction  of  the  Rayleigh  wave  arriving  from  a 
nuclear  test  in  Kazakh  from  the  coda  of  an  Aleutian  earth¬ 
quake  demonstrates  the  practical  application  of  the  technique. 
It  should  now  be  possible  to  utilize  long  period  array  data 
to  obtain  accurate  amplitude  and  phase  information  for 
small  events  which  were  previously  "hidden"  in  the  coda  of 
much  larger  events.  In  addition,  the  linear  multiple  signal 
estimator  should  make  possible  the  decomposition  of  large 


surface  waves  into  primary  and  multipath  components  on  the 


basis  of  differences  in  arrival  azimuth.  Better  estimates 
of  the  true  amplitude  and  phase  will  result  by  removal  of 
the  multipath  effects,  and  the  spectrum  and  angle  of  approach 
of  the  multipath  components  will  provide  information  as  to 
the  location  and  nature  of  the  conditions  which  give  rise 
to  multipaths. 
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Figure  1.  A  contoured  roup  of  a  2-signal  test  of  synthetic  data 
(see  page  27)  with  one  probe  point  held  fixed  while  the  other 
ranges  over  the  wavenumber  plane.  When  both  probes  occupy  the 
same  point  the  function  is  ambiguous,  its  value  varying  with  the 
direction  from  which  one  probe  approaches  the  other.  Note  the 
intersection  of  the  contours. 


AO 


Figure  2.  The  array  response  of  the  7-element  long-period 
vertical  seismic  array  at  LASA.  The  contour  interval  is 
in  3  db  steps.  The  scale  is  0.01  cycle/km  per  inch. 
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'Figure  9.  0401:53  to  0406:08  GMT,  31MAY74,  LASA  LP  array. 
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'  Figure  10.  0406:09  to  0410:24  GMT,  31MAY74,  LASA  LP  array. 
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Figure  12. 

0349: 05  to  0353:20  GMT  31MAY74,  LASA,  7  long-period 
seismometers.  The  wavenumber  spectrum  integrated  over 
the  band  0.043-0.063  Hz  after  filtering  of  secondary 
detections  derived  from  the  high-resolution  search. 

The  contours  are  given  in  3  db  intervals.  The  compon¬ 
ent  spectra  have  been  normalized. 


Figure  13. 

0353:21  to  0357:36  GMT  31MAY74,  LASA,  7  long-period 
seismometers.  The  wavenumber  spectrum  integrated  over 
the  band  0.043-0.063  Hz  after  filtering  of  secondary 
detections  derived  from  the  high-resolution  search. 
The  contours  are  given  in  3  db  intervals. 
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Figure  16. 

0406:09  to  0410:24  GMT  31MAY74,  LAS A,  7  long-period 
seismometers.  The  wavenumber  spectrum  integrated  over 
the  band  0.043-0.063  Hz  after  filtering  of  secondary 
detections  derived  from  the  high-resolution  search. 

The  contours  are  given  in  3  db  intervals. 
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